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Abstract 

We present the results of Monte Carlo simulations for the critical dynamics of the 
three-dimensional site-diluted quenched Ising model. Three different dynamics are 
considered, these correspond to the local update Metropolis scheme as well as to the 
Swendsen-Wang and Wolff cluster algorithms. The lattice sizes of L = 10 — 96 are 
analysed by a finite-size-scaling technique. The site dilution concentration p = 0.85 
was chosen to minimize the correction-to-scaling effects. We calculate numerical 
values of the dynamical critical exponents for the integrated and exponential auto- 
correlation times for energy and magnetization. As expected, cluster algorithms are 
characterized by lower values of dynamical critical exponent than the local one: also 
in the case of dilution critical slowing down is more pronounced for the Metropolis 
algorithm. However, the striking feature of our estimates is that they suggest that 
dilution leads to decrease of the dynamical critical exponent for the cluster algo- 
rithms. This phenomenon is quite opposite to the local dynamics, where dilution 
enhances critical slowing down. 
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1 Introduction 



Three-dimensional random-site Ising model (random Ising model, RIM) serves 
as a paradigm to describe influence of quenched dilution on systems which in 
a pure, undiluted, state exhibit a second order phase transition with scalar or- 
der parameter. The most prominent experimental example is given by mixed 
crystals FepZni_pF2, MnpZni_pF2 [1] however other physical realisations of 
the RIM are also possible [2]. RIM offers a unique possibility to test in the- 
ory and in simulations a change of the asymptotic critical exponents caused 
by structural randomness. Indeed, due to the Harris criterion, the new uni- 
versality class does not arise in the diluted system if the heat capacity of the 
corresponding pure system docs not diverge (i.e. if the critical exponent a < 0) 
[3]. Therefore, dilution does not change universality class of 0(m)-symmetrical 
3d systems with vector order parameter, as easy-plane or Heisenberg magnets. 
Such a universality test has been a challenge for numerous researchers in the 
on-going study since late 70-ies [1,2]. Summarized briefly, the main outcome 
of this research is that a new universality class has been found both in theory 
via renormalization group (RG) approach as well as experimentally and by 
Monte Carlo (MC) simulations [4]. 

However, the primary issue of the RIM studies performed so far concerned 
static critical behaviour and its numerical characteristics. Less is known about 
RIM critical dynamics. This paper aims to offer extensive MC simulations of 
the RIM dynamical critical properties. Moreover, as far as local and cluster 
MC algorithms correspond to different forms of dynamics, a separate task of 
our study is to compare numerical characteristics of Metropolis (local) and 
Swendsen-Wang and Wolff (cluster) dynamics for RIM. As to our knowledge 
the last question has never been addressed so far. 

The paper is organized as follows. In the next section we give a brief summary 
of available theoretical and experimental data, in section 3 we introduce the 
model and a set of obscrvables we are interested in. For the sake of complete- 
ness we briefly describe the MC algorithms in this section as well, simulation 
details are summarized in the section 4. We display and discuss the results in 
sections 5-7. 



2 Review 



In the pure (undiluted) 3d Ising model the critical slowing down, i.e. an in- 
crease of the relaxation time r as the critical point is approached, is gov- 
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erned by the universal dynamical critical exponent z: 

T~|r-r,|— , (1) 



with the correlation length critical exponent u. For isotropic systems, the 
dynamical exponent is related to the pair correlation function critical exponent 
r] via z — 2 + crj [11], where c is a (d-dependent) constant. The numerical 
value of the exponent i] being small, the value of the dynamical exponent z 
for the 3d Ising model slightly differs from 2. Typical numbers are z = 2.1(1) 
(experiment, FeFs, [5]), z = 2.032(4) (MC, [6]), z = 2.017 [7], z = 2.012 [8] 
(RG theories). 

Below we briefly review experimental, theoretical and MC studies performed 

so far to analyse how the relation (1) holds for the RIM and, in particular, 
to answer the question how the dynamical critical exponent z is influenced by 
dilution. Numerical values of z that follow from our review are collected in 
Table 1. 

Experiments. There are only three independent studies of the RIM critical 
dynamics we are aware of. All three concern an antiferromagnetic uniaxial 
crystal FeF2 diluted by its non-magnetic isomorph ZnF2. The resulting sub- 
stance FepZni_pF2 was analysed by two different techniques: Mossbauer spec- 
troscopy (studying the dynamical hne broadening of the Mossbauer spectra) 
[9,10] and by the spin-echo neutron scattering (analyzing the time dependent 
spin correlation function) [5]. The earlier studies [9,5] lead to conclusion that 
dilution causes a decrease of z: it was found to vary from 2.1 to 1.5 when mag- 
netic atom concentration p varied from 1 to 0.46 [9], whereas Ref. [5] reports 
z — 1.7(2) for Feo.46Zno.54F2 and explains it in the frames of conventional van 
Hove theory: z = 2 — rj. However, a later experiment [10] brings about the 
value z = 2.18(10) for Fco.9Zno.1F2 being in a good agreement with available 
RG results, as we will see below. 

Theory. Theoretical analysis of the RIM critical dynamics is due to the RG ap- 
proach. The majority of work was devoted to analysis of the pure relaxational 
dynamics with non-conserved order parameter described by the Langevin 
equation of motion, model A dynamics in the classification of Ref. [11]. Be- 
ing realized via MC simulations, it corresponds to the single spin Metropolis 
dynamics [12]. Change of the RIM equations of motion by couphng to the 
diffusive dynamics of a conserved scalar density (energy density) does not 
change the asymptotic critical behavior [13]: model C and model A belong to 
the same dynamical universality class for the RIM. However, due to crossover 
effects, the effective critical behaviour essentially differs for these two forms of 
dynamics [14]. 
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Table 1 

The dynamical critical exponent of RIM z as defined in experiments, theory and 
MC simulations. The methods are given in the following notations, experiments: MS, 
Mossbauer spectroscopy; SENS, spin-echo neutron scattering; theory: e^/^, first non- 
trivial order e-expansion; MRG, massive RG at (i = 3; MS, minimal subtraction RG 
at d = 3; Metropolis MC simulations: FSS, finite-size-scaling; DRG, dynamical RG; 
OE, out-of-equlibrium short-time dynamics. See the text for a whole description. 



Reference 


Method 


Peculiarities 


z 


Barrett et al., 1986 [9] 


MS 


FepZni_9F2, 








0.46 < p < 1 


1.5<z< 2.1 


Belanger et al., 1988 [5] 


SENS 


Feo.46Zno.54F2 


1.7(2) 


Rosov et al., 1992 [10] 


MS 


Feo.9Zno.1F2 


2.18(10) 


Grinstein et al., 1977 [15] 


£1/2 




2.336 


Prudnikov et al., 1992 [17] 


MRG 


2 loops 


2.237 


Janssen et al., 1995 [20] 


MS 


2-3 loops 


2.18 


Prudnikov et al., 1998 [18] 


MRG 


3 loops 


2.165 


Blavats'ka et al., 2005 [8] 


MS 


2 loops 


2.172 


Prudnikov, 1992 [32] 


DRG 


L = 48, ZM 


z{p = 0.95) = 2.19(7) 








z{p = 0.8) = 2.20(8) 


Heuer, 1993 [31] 


FSS 


L = 60, 2;|M|,int 


2.4(1) 


Parisi et al., 1999 [34] 


OE 


L = 100, z^ 


2.62(7) 


Schehr et al., 2005 [26] 


OE 


L = 100, ZM 


2.6(1) 



The pioneering work [15] analysed by £ = 4 — d expansion the critical dynamics 
of the m-vector model with quenched random impurities and non-conserved 
order parameter. For the RIM case, m = 1, new dynamical universality class 
was found with the two loop value of the exponent z = 2 + ^6e/53. Being 
analysed at 0? = 3 naively by a simple substitution e = 1 this yields z ~ 
2.336 and essentially differs from z of the pure 3d Ising model quoted at 
the beginning of this section. However, the RG expansions are known to be 
asymptotic at best and resummation is needed to get reliable data on their 
basis [16]. Further results were obtained by the massive RG approach directly 
at (i = 3 with subsequent resummation of resulting expansions: the two-loop 
value of the exponent reads z — 2.237 [17], the three-loop one is 2; = 2.165 
[18]. Currently, due to essential technical difficulties, dynamical RG functions 
of RIM are known only within two-loop accuracy in the minimal subtraction 
RG scheme and within three loops in the massive RG approach at d = 3. In 
the minimal subtraction renormalization, the static RG functions were taken 
in three loops and combined with the two-loop expansions for the dynamic 
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ones. These gave the following e-expansion for the dynamical exponent: z — 
2 + 0.336v/i(l - 0.932v^) with the naive estimate z = 2.023 [19]. Being 
improved by the resummation of the static RG functions, the estimate reads 
z — 2.18 [20]. The last value is close to the other estimate, obtained from the 
resummation of the two-loop minimal subtraction RG functions directly at 
d = 3: z = 2.172 [8]. 

Theoretical studies mentioned above concerned dynamic criticality associated 
with equilibrium fluctuations. Recently, it was shown [21] that non-equilibrium 
relaxation at short times possess scaling features as well. In particular, if a 
system is suddenly quenched from high temperatures to the critical one and 
then released to the dynamic evolution of model A, the evolution at short 
times is governed by scaling laws. Both dynamical and static critical exponents 
can be extracted from the scaling. Short-time dynamics of the RIM at non- 
equilibrium critical relaxation was analysed if Refs. [25,19,26,27]. 

Let us mention a related problem, where an influence of quenched disorder 
on RIM critical dynamics was examined theoretically. These are studies of an 
effect of extended impurities on RIM critical dynamics [13,28,14]. As far as 
the presence of extended (long-range correlated) impurities changes the static 
universality class of RIM, the dynamical critical behaviour is found to differ 
from those of the RIM with point-like uncorrelated disorder. 

MC simulations. Essential progress in MC simulations of static critical phe- 
nomena is due to the application of cluster algorithms [29,30]. In particular, 
they allowed to obtain precise values of the RIM static critical exponents [2]. 
It is to be emphasized here, that whereas the cluster algorithms were spe- 
cially designed to lead to the same static critical behaviour as the single-spin 
Metropolis algorithm [12], it is not the case for dynamics. It is the Metropolis 
algorithm (due presumably to its locality), which leads to the same value of 
dynamic critical exponent as the one observed for RIM experimentally in Refs. 
[5,9,10] and analysed theoretically in Refs. [15,17,18,19,20,8]. In cluster algo- 
rithms, as it follows already from their name, the whole clusters of spins are 
flipped, which gives origin to the non-local dynamics. In its turn, the last is 
characterized by its own dynamical critical exponents [22,23,24]. As far as the 
cluster algorithms were introduced to overcome the critical slowing down, the 
corresponding autocorrelation times are characterized by weaker singularities 
as those of local dynamics: ^iduster < 2;iocai- 

Let us note, that the theoretical RG calculations assume a single dynamical 
critical exponent z for the relaxation times of different observables. In the MC 
simulations one typically flnds, that the autocorrelation time of different ob- 
servables is characterized by different (effective) exponents, which are expected 
to coincide in the asymptotics. Therefore, when we give the MC values of z 
in Table 1 we specify also the physical observable for which it has been mea- 
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sured. Already the first MC study of tlie RIM single-spin critical dynamics 
revealed dynamical scaling behaviour with a concentration-dependent criti- 
cal exponent z [31]. For small dilution the following numbers were reported: 
z{v = 0.95) = 2.15(1); z{v = 0.9) = 2.23(1); z{p = 0.8) = 2.39(1). The concen- 
tration dependence of z was explained by crossover. In the region of concentra- 
tions p ~ 0.8 the slope of the crossover function was found to change its sign, 
therefore the correction-to-scaling terms were minimal. This allowed to arrive 
to the conclusion about an asymptotic value of the dynamic critical exponent 
z — 2.4(1) [31]. Independently, Metropolis dynamics of the RIM was analysed 
in Ref. [32] by a combination of MC and dynamical RG [33] techniques. Again 
the concentration-dependent exponents were found with a different conclusion, 
however: a hypothesis of RIM step-like universality was proposed. According 
to the hypothesis, the asymptotic critical exponents remain unchanged only 
within certain concentration region. For a small dilution the values of z prac- 
tically did not differ: z{p = 0.95) = 2.19(7); z{p = 0.8) = 2.20(8) [32] and are 
compatible with those obtained by a finite-size-scaling technique in Ref. [31]. 

Other estimates come from MC simulations of the out-of-equilibrium RIM dy- 
namics. Here, taking into account correction-to-scaling, value z — 2.62(7) was 
extracted from the time dependence of the out-of-equilibrium susceptibility 
x{t) (with the leading dynamical correction-to-scaling exponent u — 0.50(13)) 
[34]. This value was further supported by the out-of-equilibrium simulations 
of Ref. [26], where the value z = 2.6(1) was extracted from the scaling of the 
spin-spin autocorrelation function. 

As it was noted above, the MC cluster algorithms provide different type of 
dynamics and therefore their scaling exponents can not be compared straight- 
away with those of single-spin local dynamics summarized in Table 1. More- 
over, currently there is no field theory available to predict the critical exponent 
value for cluster dynamics even for the pure (undiluted) spin models. Study 
of such dynamics constitutes a separate task and certain analytical and nu- 
merical work has already been done for the pure models [22,23,24]. As to our 
knowledge, no results for the RIM have been obtained so far. An exception 
is Ref. [35], where an effective (concentration dependent) critical exponent z 
was obtained for the Swendsen-Wang cluster algorithm for 3d random-bond 
Ising model. For a small bond dilution an estimate reads: z{p — 0.7) = 0.41 
[35]. 



3 Observables and MC algorithms 

In our paper we consider the 3d Ising model with non-magnetic impurities 
randomly distributed over the system. The Hamiltonian of this model on the 
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cubic lattice has the following form 

H ^ -J^CiCjSiSj, (2) 

{ij) 



where {ij) denotes the summation over the nearest neighbour sites of the 
lattice, Cj = 1 if the i-th site is occupied by a spin and Cj = otherwise, the 
Ising spins Si take on the values +1 or —1. The spins interact via an exchange 
coupling J, which is positive. Occupied sites (q = 1) are considered to be 
uncorrclatcd, randomly distributed and quenched in a fixed configuration. For 
every observable discussed below, first the Boltzmann average with respect 
to the spin subsystem is performed for the fixed disorder realisation, then 
the averaging over different disorder realisations is performed. We will use 
the following notations: Boltzmann average over the spin subsystem will be 
denoted by angular brackets ((. . .)) whereas the over bar (. . .) will stand for 
the averaging over the disorder realisations. The number of all sites is = 
and the number of sites carrying a spin is Np. The concentration of spins is 
defined therefore as p = Np/N. 



3. 1 The properties of interest 



For a given disorder realisation, an average value of an observable (C) at 
temperature T can be computed in the canonical ensemble from its values O 
for given spin configurations: 

(e») = |spOe-^«, (3) 



where (3 — (kT) ^, Z is the partition function 

Z = Spe-'^^, (4) 



and trace in (4) is taken over the spin degrees of freedom. In the course of the 
MC simulation each spin configuration is generated with its proper Boltzmann 
weight already (more detailed description of the algorithms is given in the 
next section) , hence the thermodynamic average is the simple average over all 
generated configurations. If A^steps is the total number of productive MC steps 
used for the averaging, then 

= E (5) 

-'V steps conf 
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Sum in (5) spans all spin configurations in which the (spin configuration de- 
pendent) observable O is measured. 



In an ideal case of uncorrelated, statistically independent configurations, the 
total error in defining (O) can be evaluated as 



^ -^steps ( -^steps 1) 



where 60 = O — (O). In practice, however, the correlation between differ- 
ent spin configurations exists as the result of particular MC scheme. This 
correction can be characterized via the (disorder dependent) autocorrelation 
function [36]: 

_ {60{to + 6t)60it,)) 
^""^^'^ - {50{t, + 6t)){60{t,)y 



where to is some time origin. 

At times large enough Co{St) decays exponentially according to the Debye 
law 

St /to,, 



where To,exp is the exponential autocorrelation time, obtained for quantity 
O and a is a constant. Time To,exp defines a time scale at which the con- 
figurations generated in a course of the MC simulation can be assumed as 
uncorrelated. Hence, in the simulation run of length A'^stcps MC steps, only 
„"^°*''''° configurations are considered to be statistically independent. 



Besides the r^^exp, the relaxation of an observable O is characterized by the 
integrated autocorrelation time To,int- It is defined via 

-1 oo 

roM-7;+T.Co{St)). (9) 

^ St=l 



In practice, To,mt is evaluated by introducing a maximum cutoff in the sum (9) 
and it may be shown that both autocorrelation times coincide To,int — T"C',exp 
only in the limit when this cutoff goes to 5t — > oo, otherwise To int < ^Oexp 
[38]. 

Let us specify now the observables we will be interested in during MC simula- 
tions. In this study we concentrate on the internal energy S, the magnetization 
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M. and absolute value of the magnetization \M.\ per site, defined as 

£ = - J-^ ^ CiCjSiSj, (10) 

P (ij) 

M=^T.C^S^, (11) 
P i 

1-^1 = 1^1 E^^-^^l- (12) 



Prom these observables {O) we compute the following expectation values (O): 



E^{£), M^{M), |M| = (|M|). (13) 



The divergency of the autocorrelation time for any of these quantities as the 
Tc is approached, Eq. (1), for a finite system of size L is reflected in a power 
law scaling of r with L. Taken that both exponential and integrated auto- 
correlation times are defined during simulation, their scaling is governed by 
corresponding exponents: 



ro,exp ~ L'o,"^^, (14) 

T0,int - (15) 

with O being any of the expectation values E, M, \M\ computed in simula- 
tions. 

Note, that the theoretical RG calculations assume a unique dynamical critical 
exponent z for the relaxation times of all observable quantities (cf. theoret- 
ical estimates for z in Table 1). This might not hold for the scaling of the 
experimentally defined autocorrelation times (14), (15). For the pure system, 
it is believed that the scaling of autocorrelation times for the energy-like ob- 
servables is described with the same dynamical critical exponent [37], the last 
might differ from that for the susceptibihty-like observables [24] . 

As it was noted already above, local and cluster MC algorithms give rise to 
different forms of critical dynamics. Therefore, they are described by different 
scahng exponents (14), (15). Indeed, the large clusters of equally oriented 
spins are started to be formed in the vicinity of a critical point. Therefore, 
most of the spin update attempts made by any of local algorithms (e.g. by 
the Metropolis one) are wasted, and, as result, the generated configurations 
are highly correlated. Both the correlation time r and the critical index z are 
large and system moves in a configurational space inefficiently. This poses the 
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severe difficulty for the calculation of the static critical exponents, but, in fact, 
reflects the real dynamics in the system near the critical point. 

The non-local, cluster algorithms introduced by Swendsen-Wang and Wolff 
consider cluster pseudo-dynamics of the system by attempting to flip the whole 
cluster(s) at once. The primary goal is to reduce the effective autocorrelation 
time and, therefore, to improve the statistical sampling of generated config- 
urations. At the same time, these algorithms bring the system into different 
dynamical class and an analysis of the different dynamics governed by dif- 
ferent MC algorithms is the main goal of this study. Below, for the sake of 
completeness, we briefly describe main steps of the MC algorithms under con- 
sideration. 

3.2 MC algorithms 

The Metropohs algorithm is the simplest and historically the first MC algo- 
rithm [12]. It utilises the preferential sampling of the configurational space, 
where each spin configuration is generated with appropriate Boltzmann weight. 
The technical difficulty of generating all the new configurations independently 
is overcame by using instead the recipe how to produce each new configura- 
tion from the previous one. Spins are randomly selected and flipped with a 
probability P = min{l, e~^^^) where A/J is the energy difference between 
the old and the new configurations. A Markovian chain of configurations is 
thus produced and defines the pseudo dynamics of the system. 

Locality of this algorithm is a serious drawback near the critical point, where 

large correlated clusters of spins are emerging. As the result, the vast amount 
of single spin flips are rejected, therefore the conflgurations generated are 
highly correlated and the system moves in a phase space inefficiently. 

Cluster algorithms were designed to overcome these difficulties. They are based 
on the idcntiflcation of clusters of sites using a bond percolation process con- 
nected to the spin conflguration of the magnetic system. All spins of the clus- 
ters are then independently flipped. 

In the case of the Ising model (or more generally the Potts model), the percola- 
tion process involved is obtained through the mapping onto the random graph 
model, flrst addressed by Fortuin and Kasteleyn [40]. In the Swendsen-Wang 
algorithm [29], a cluster update sweep consists of the following steps: depend- 
ing on the nearest neighbour exchange interactions and site occupations, assign 
value to a bond between sites i and j with probability Pij = 1 — ^-'^PJciCj ^ then 
identify clusters of spins connected by active bonds, and eventually assign a 
random value to all the spins in a given cluster. The spin system at criti- 
cality is mapped into a bond percolation problem at the percolation thresh- 
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old. It results in the producing of clusters of arbitrary large sizes, therefore 
the Swendsen-Wang algorithm samples the configurations in a critical region 
much more efficiently. This is reflected in its rather small value of the dynami- 
cal critical exponent z. Some disadvantage is an extra computer time required 
to split-up the system into clusters. 

Wolff introduced a single cluster algorithm [30] which otherwise is much sim- 
ilar to the Swendsen-Wang one. A spin is randomly chosen, then the cluster 
connected with this spin is constructed and all the spins in the cluster are 
updated. Note that in this scheme, the flip of one cluster updates only the 
spins belonging to this cluster and therefore produces only a partial update of 
the system. To match the same time scale as in the Metropolis and Swendsen- 
Wang algorithms, the time scale of the Wolff algorithm should be corrected 
by a factor c = iciuster/-t'^, where ^cluster is the average size of flipped clusters. 



4 Simulation details 

In the rest of the paper we study RIM critical dynamics governed by the three 
different MC algorithms described in the previous section. The MC simulations 
are performed for a range of system sizes up to L = 96 with periodic boundary 
conditions. The concentration of magnetic sites was flxed at p = 0.85. This 
choice is based on the previous findings that the correction-to-scaling terms 
are minimal at concentrations p ~ 0.8 [31,41]. Therefore, we do not account 
for these terms. The simulations are performed at the temperature that corre- 
sponds to the critical temperature of the infinite system and was taken equal 
to ficJ = 0.2661922(83) according to our previous findings [42]. 

Table 2 



CPU time (in seconds) used for performing 1000 MC steps and 10 disorder realisa- 
tions, 



L 


10 


12 


16 


24 


32 


48 


64 


96 


Metropolis 


3.33 


5.81 


13.77 


50.17 


119.26 


403.69 


957.72 


3207.54 


Sw.-Wang 


3.60 


6.18 


14.60 


55.10 


130.78 


442.64 


1092.76 


3682.48 


Wolff 


1.25 


2.05 


4.53 


18.11 


43.39 


139.81 


298.33 


970.05 



Typically, we averaged over A'^dis = 10^ disorder realisations for each lattice 
size. All runs were started from a random configuration of empty and occupied 
spin sites. At first, we run 250rE,int MC sweeps for thermal equilibration and 
then the production run of lO^rg^int MC sweeps was conducted. For all cases 
that was quite sufficient for the accurate description of the long-time behaviour 
of the autocorrelation functions. 

As a basis for the random number generator we take minimal random num- 
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ber generator Rani of Park and Miller with Bays-Durham shuffle and added 
safeguards, described in Ref. [43]. 

We used the workstations cluster of the ICMP based on Athlon MP 2200+ 
processors. The typical simulation time per 1000 MC sweeps and 10 disorder 
reahsations with data records are shown in the Table. 2. 



5 Autocorrelation times 



When performing the MC simulations on a system of N spins, the time unit 
can be related to the number of MC sweeps (MCS), where during one sweep 
on average N spins are updated. This convention is also used in our study. 
A quantitative analysis of the autocorrelation times involves an evaluation of 
the autocorrelation function for various properties of the system [36] . As was 
mentioned, in this study we concentrate on the autocorrelations of the energy 
S, the magnetization Ai and the absolute value of magnetization \ Defined 
in (7), the autocorrelation function can be rewritten as: 



{50{to + 6t)){60{to)) 

{Ojto + 5t)0{to)) - {Ojto + 6t)){0{to)) 

^mto+sm - {o{to+sm{{o{tor) - mom 



(16) 



In the thermodynamic limit {O{to + 6t)) is equal to {O{to)), and one arrives 
to the more usual expression 

^ .... ^ mo)0{to + 6t))- {Ojto)) {Ojto)) 

' {0{to)0{to))-{0{to)){0{to)) ' ^ ^ 



where 0{t) is the instant value for the property of interest at certain time t, 
(to is some time origin, 5t is the time elapsed since the time origin to). The 
averaging over a large number of time origins to is needed to smoothen up the 
Co{5t) at large 5t. As an example, a typical behaviour of the energy-energy 
autocorrelation functions for given disorder realisations is shown in Figs. 1-3. 

Table 3 

The average cluster size /cluster and factor c for Wolff MC algorithm 



L 


10 


12 


16 


24 


32 


48 


64 


96 


^cluster 


149 


216 


385 


856 


1526 


3393 


6241 


15034 


C 


0.175 


0.147 


0.111 


0.073 


0.055 


0.036 


0.028 


0.020 



12 



C^(St) 1 



0,1 



0,01 



1E-3 



Metropolis algorithm 




500 1000 1500 

5t 



Figure 1. The log-linear plot for the energy-energy autocorrelation function Cs{6t), 
the Metropolis algorithm, L = 64. Bold line: measured value. Dashed line: fit to the 
exponential decay (8) with the autocorrelation time T^^exp = 507.8. 



q(st) 



0,01- 



1E-3- 



Swendsen-Wang algorithm 



20 



40 



60 



St 



Figure 2. The log-linear plot for the energy-energy autocorrelation function Cs{6t), 
the Swendsen-Wang algorithm, L = 64. Bold line: measured value. Dashed line: fit 



to the exponential decay (8) with the autocorrelation time T£ 



exp 



9.04. 



As was already mentioned above, the time scale of the Wolff algorithm should 
be accounted for the average cluster size to be compared correctly with the 
dynamics of other algorithms. To this end, for each lattice size and for a 
given disorder realisation we calculated the size of flipped cluster, skipping 
the first 250rs^int MC steps for thermal equilibration. Then we performed 
configurational averaging of the updated cluster size. The scaling factor c 
was calculated as c = /ciustcr/A'p, where /cluster is the cluster size averaged 
over different disorder realisations. In Table 3, the average cluster size and 
the scaling factor c are presented. One could explain the behaviour of the 
scaling factor c by the following considerations. For the smaller system sizes, 
the simulation temperature (which is equal to T^) is much lower than the 
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Figure 3. The log- linear plot for the the energy-energy autocorrelation function 
C£{5t), the Wolff algorithm, L = 64. The measured value and the fit to the expo- 
nential decay (8) with the autocorrelation time r^-^exp = 2.63 are indistinguishable 
within the scale chosen. 

effective critical temperatures T^, the system is at T < so that the typical 
cluster occupy larger part of the system. With the increase of the system size, 
the simulation temperature is getting closer to the and the average clusters 
being flipped are decreasing in size. 

To calculate the integrated autocorrelation time To,int the expression (9) is 
used. One should note that the error for the autocorrelation function is always 
larger at long times, where the data are averaged over less intervals. The 
compromised accuracy for the integrated autocorrelation time can be achieved 
by using certain time cutoff Stmax, typically of the order of Stmax > 6rj„( [38]: 

roMStmax) = ^ + E CoiSt)), (18) 

^ St=l 



The trailing part of the autocorrelation function for 6t > 6tmax can be ap- 
proximated by an exponential function. As the result, the flnal expression 
reads 



1 Stmax OO 

TOM = 7;+ i: Com) + « E e-^*/*°.- = 

5t=l 5t=5tmax+l 

= roM^tmax) + a- e-^t^a./ro,exp_ (19) 



In this study, we employed the following scheme. For each disorder realisa- 
tion the autocorrelation function CoiSt) has been evaluated flrst. To calculate 
first term in (19) we use (18) with condition of the cutoff: Stmax > Sro^nt- 
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To evaluate the second term, one has to estimate To,exp- Whereas the pure 
systems display asymptotic behaviour dominated by a single relaxation time, 
the distinct feature of the autocorrelation function of disordered systems is 
that there is a whole spectrum of autocorrelation times in the crossover region 
and it is reflected in the curvature of the autocorrelation functions in the log- 
log plot. This feature has been noted already in Ref. [31]. Therefore, we plot 
the Co{St) in log- log scale and estimated To,cxp and a from the straight line 
region in a window To,mt{Stmax) to STo^mti^tmax)- The averaged over disorder 
realizations exponential autocorrelation time ro,exp is given in Table 4 . 

Table 4 

Exponential autocorrelation times ro,exp of RIM at the critical temperature of in- 
finite system for different lattice sizes L measured in MC sweeps for Metropolis, 
Swendsen-Wang and Wolff MC algorithms. 



L 


Metropolis 




Swendsen-Wang Wolff 






TE,exp 




''"|M|,exp 


TM,exp TE,exp ''"|M|,exp TE,exp 


''"|M|,exp 


10 


8.94(63) 




9.34(88) 


6.80(1.80)40 3.76(19) 3.72(20) 1.53(17) 


1.23(10) 


12 


1.33(11). 


10 


1.40(16)' 


10 1.02(28)402 4.16(23)4.12(21)1.64(17) 


1.30(9) 


16 


2.49(24) ■ 


10 


2.69(44)' 


10 1.98(63)402 4.82(27)4.80(36)1.80(21) 


1.41(10) 


24 


5.99(37)' 


10 


6.51(95)' 


10 5.10(19)402 5.87(38)5.92(37) 2.06(24) 


1.54(10) 


32 


1. 12(08)' 


102 


1.24(28)' 


102 1.02(57)40^ 6.78(52) 6.71(42) 2.20(28) 


1.65(11) 


48 


2.71(25)' 


102 


2.86(40)' 


102 2.40(1.30)40^ 8.15(69) 8.02(52) 2.51(30) 


1.75(10) 


64 


4.65(48)' 


102 6.31(97)' 


102 4.78(2.25)40^ 9.19(85) 9.03(76) 2.35(21) 


1.88(11) 


96 


1.22(14)' 


102 


1.44(26)' 


10=^1.20(0.52)40^ 10.7(9) 10.5(9) 2.79(35) 


1.99(16) 



In order to calculate the error bars for the exponential and integrated relax- 
ation times we use the blocking method. We divide all set of autocorrelation 

times (each corresponding to a separate replica) into n blocks so that each 
block contains values of the autocorrelation times. We obtain the average 
value and standard error due to formulas (5) and (6). Then we do simple av- 
eraging over n blocks. We do not give results for TM,exp and TM,int for cluster 
methods, because the correlation of M are absent. 



6 Critical exponents 

Having found autocorrelation times for different observables as functions of lat- 
tice size one can extract via Eqs. (14), (15) the values of dynamical exponents 
for each of the algorithms considered. Let us start with the Metropolis algo- 
rithm. Log-log plots for the integrated and exponential autocorrelation times 
for E, M, and \M\ are shown in Fig. 4. We used a hnear square interpolation 
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Table 5 

Integrated autocorrelation times To,mt of RIM at the critical temperature of infinite 
system for different lattice sizes L measured in MC sweeps for Metropolis, Swendsen- 
Wang and Wolff MC algorithms. 



L Metropolis 




Swendsen-Wang 


: Wolff 


TE, int T\M \,int 


T M,mt 


TE,int ''"|M|,int 





10 6.14(30) 9.22(49) 6.6(1.6)-10 3.57(16)3.22(17) 1.33(8) 0.95(4) 

12 8.47(49) 1.36(8)40 9.8(2.5)-10 3.92(18)3.50(19) 1.42(9) 0.98(4) 

16 1.44(10)40 2.53(20)40 1.87(53)40^ 4.47(23)3.95(25) 1.55(12) 1.00(4) 

24 3.10(20)40 6.13(48)40 4.74(1.50)402 5.33(31)4.68(27) 1.75(15) 1.02(4) 

32 5.51(40)40 1.16(14)40^ 9.05(3.41)40^ 6.05(36) 5.14(34) 1.87(18) 1.04(4) 

48 1.22(12)40^ 2.71(28)402 2.17(0.82)40^ 7.11(47) 6.07(43) 2.08(20) 1.04(4) 

64 2.15(20)40^ 5.59(62)402 3.80(1.44)40^ 7.81(51) 6.71(50) 2.16(18) 1.07(3) 

96 4.84(51)402 1.24(15)40^ 9.57(3.45)40^ 9.21(61) 7.65(59) 2.33(25) 1.10(4) 

to extrapolate the t{L) dependencies to get the values of the exponents. To 
extrapolate data obtained from the Metropolis algorithm all eight data points 
were used, whereas for Swendsen-Wang and Wolff algorithms only five last 
data points (the largest system sizes) were considered. 

Log-log plots for the autocorrelation time for cluster algorithms are shown 
in Fig. 5 (Swendsen-Wang algorithm) and Fig. 6 (Wolff algorithm). As noted 
above, not all autocorrelation times are well-defined for the cluster algorithms. 
In particular, for the Swendsen-Wang algorithm we were able to define auto- 
correlations for E and \ M\ (and not for the magnetization per site M), whereas 
for the Wolff algorithm only r^'^int, te,cxp, and T|m|,cxp were well-defined. How- 
ever, to estimate value of the exponent ZE,exp we have to discard data for 
L — 64. As one can see from the plot in Fig. 6b an appropriate data point 
when included does not lead to a reasonable linear approximation. 

Numerical values of the exponents are given in the first line of Table 6. With 
the exception of the exponent for the integrated energy autocorrelation time 
TE,int, the rest of the exponents definitely are close to the value 2; ~ 2.2. Com- 
paring this value with the data of Table 1 one sees, that it is in a reasonable 
agreement with the theoretical estimates of [15,17,18,20,8] as well as with the 
experimental result of [10] and data of MC simulations [31,32]. Note how- 
ever, that estimates from the out-of-equilibrium MC simulations give rather 
different value z ~ 2.6 [34,26]. 

A rather striking feature of the dynamical exponents of cluster algorithms is 
that dilution leads to decrease of the exponents, as compared to the pure 3d 
Ising model. Indeed, the value for the integrated dynamical critical exponent 



16 




2,0 2,5 3,0 3,5 4,0 4,5 ln(L) 2,0 2,5 3,0 3,5 4,0 4,5 ln(Z) 




(e) (f) 

Figure 4. Integrated (left column) and exponential (right column) autocorrelation 
times as functions of L for the Metropolis algorithm. a,b: energy autocorrelation, 
'TE;mti T^.exp; Cjd: magnetization autocorrelation, tm, iat-, TM,exp] e,f: absolute value 
of magnetization autocorrelation, r|A/|_int, |^exp- 

of the Swendsen-Wang algorithm for the "energy like" observables recently 
calculated in Ref. [24] reads: ^^int = 0.459(30). In the same study, the dy- 
namical critical exponent associated to the exponential autocorrelation time 
was found to be Zexp = 0.481. Both values exceed those found by us for the 
RIM, see the second line of Table 6. Similar tendency to the dilution induced 
decrease of the Swendsen-Wang dynamical critical exponent was observed re- 
cently for the random-bond 3d Ising model [35]. There, the value z = 0.41 
for the bond concentration p = 0.7 was found, which again is smaller than its 
counterpart for the pure 3d Ising model [24]. 
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^ e.exp ' 

2,4 - 



2,2 - 

2,0 - 

1,8 - 

1,6 - 

1,4 - 
1,2 



z =0,44(7)- 



t 



2,0 2,5 3,0 3,5 4,0 4,5 ln(L) 

(b) 



Vi.exp ) 




Figure 5. Integrated (left column) and exponential (right column) autocorrelation 
times as functions of L for the Swendsen-Wang algorithm. a,b: energy autocorre- 
lation, r£;_int, T£;,exp; c,d: absolute value of magnetization autocorrelation, T|jvf|,int) 

''"IM|,exp- 

Table 6 

RIM dynamical critical exponents for different MC algorithms, obtained from data 
at the critical temperature of the infinite system. Tabs. 4-5 







^i?,exp 


^|Af|,int 


^|Af|,cxp 


ZM,mt ZM,cxp 


Metropolis 


1.92(4) 


2.16(4) 


2.18(4) 


2.23(6) 


2.21(16) 2.29(19) 


Swendsen-Wang 


0.39(6) 


0.44(7) 


0.36(6) 


0.42(7) 




Wolff 


0.21(8) 


0.22(12) 




0.19(6) 





We re-analysed simulation data, considering them at the critical temperature 
of a finite system of size L, Tc{L). The latter may be calculated in different 
ways, being defined by the maximum of different observables. In Table 7, we 
give the values of RIM dynamical critical exponents at Tc{L) obtained from 
the maximum of magnetic susceptibility. The critical temperature was taken 
from our previous study [42]. As one can see comparing Tables 6 and 7, the 
crossover effects do not influence data essentially. 

One more question worth discussing is whether the relations between cluster 
algorithms dynamical exponents and the static exponents observed for the 
pure systems [23,22] hold for the diluted ones. Indeed, for the pure Ising model. 
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^ E.mt ■' 



0,8 



0,6 - 



0,4 



0,2 



Z . =0,21(8) 



ln(^, ) 

^ e,exp ' 




z =0,22(1 2y 



2,0 2,5 3,0 3.5 4,0 4,5 ln(i) 

(a) 



2,0 2,5 3,0 3,5 4,0 4,5 ln(L) 

(b) 



0,6 



0,4 




z,„, =0,19(6) 



2,0 2,5 3,0 3,5 4,0 4,5 ln(L) 

(c) 

Figure 6. Autocorrelation times as functions of L for the Wolff algorithm. a,b: 
energy autocorrelation, r£;,int, TE,exp\ c: absolute value of magnetization autocorre- 
lation, r|A/|^exp- 

Table 7 

RIM dynamical critical exponents for different MC algorithms, calculated at the 
critical temperature of the finite size system Tc{L). 



ZE;mt ZE,cxp ^|Af|,int Z\M\,cxp ZM,int 



ZM. 



exp 



Metropolis 

Swendsen-Wang 

Wolff 



1.99(3) 2.22(3) 2.19(5) 2.22(7) 2.18(12) 2.23(16) 
0.35(5) 0.39(6) 0.32(6) 0.40(6) - 
0.16(8) 0.16(9) - 0.14(5) - 



the Coddington-Ballie conjecture holds [23], stating that the Swendsen-Wang 
dynamical critical exponent ^f^nt is defined via static critical exponents for 
magnetization and correlation length: 



-EAnt 



(20) 



It is worth here to note, that whereas the static critical exponents for RIM 
numerically differ from those of the pure 3d Ising model (cf. theoretical RG 
estimates (3 = 0.349(5) and u = 0.678(10) [45] for RIM with (3 = 0.3258(14) 
and u = 0.6304(13) [46] for pure 3d Ising model) their relation remains al- 
most unchanged. For the numbers given above, l3/v = 0.515(15) for RIM and 
13/ V = 0.517(3) for the 3d Ising model. Therefore, a change in the value of the 
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Swendsen-Wang dynamic critical exponent upon dilution serves as an evidence 
that the relation (20) does not hold for RIM. 

Another empirical relation found in Ref. [23] for the pure Ising model connects 
the dynamical critical exponent of the Wolff algorithm z^^^^t '^^^^ the static 
ones: 

^^int = «/^- (21) 



As far as a < /3 for the 3d Ising model, comparison of Eqs. (20) and (21) leads 
to the inequality: 



Eq. (21) does not hold for the diluted systems (where the heat capacity crit- 
ical exponent is negative). However, the inequahty (22) still holds, as one 
can see, comparing data of Table 6 for the Wolff and Swendsen-Wang algo- 
rithms. Again, the value of dynamical critical exponent for the diluted system 
is smaller than its counterpart for the pure one. Wolff algorithm dynamical 
critical exponent of the 3d Ising model found in different simulations read: 
zem = 0.28(2) [47]; 0.44(10) [48]; 0.33(1) [23], all numbers exceeding those 
given for the Wolff case of RIM in Table 6. 



7 Conclusions and outlook 



In this paper, we have studied dynamical critical behaviour of the 3d random- 
site Ising model (RIM) originated from different MC algorithms. We consid- 
ered the local single-spin Metropolis algorithms as well as Swcndscn-Wang 
and Wolff cluster algorithms. Giving origin to an equivalent static critical be- 
haviour, all three algorithms correspond to different forms of dynamics. A 
comparison of numerical characteristics of Metropohs (local) and Swendsen- 
Wang and Wolff (cluster) dynamics for RIM was achieved by calculation of the 
integrated and exponential autocorrelation times for RIM energy and magne- 
tization. 

The local update Metropolis algorithm corresponds to the pure relaxational 
single-spin dynamics with non-conserved order parameter and finds its theoret- 
ical description as the model A critical dynamics [11]. There exist RG analysis 
of critical dynamics for the RIM with such type of relaxation [15,17,18,20,19,8]. 
It assumes, however, a single dynamical critical exponent z for the relaxation 
times of different observables. Although the perturbation theory series are 
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divergent and only the lowest non-trivial order calculations have been per- 
formed so far, being resummed appropriately, all available theoretical data 
are coherent with an estimate z ~ 2.2(1) (see Table 1 of our paper). This is 
further supported by the latest experimental observation we are aware about, 
z — 2.18(10) [10]. Whereas initial MC simulations gave estimates of z along 
with the above value [31,32], recent simulations of Refs. [34,26] favour an esti- 
mate z ~ 2.6(1). Our results for the values of RIM dynamical critical exponent 
for a local dynamics are summarized in the first line of Table 6. The essential 
features of the discussion remain unchanged, when one considers calculations 
at the critical temperature of the finite size system, Tc(L), Table 7. Except of 
the exponent for the energy integrated autocorrelation time, ZE,int, our data 
supports an estimate z ~ 2.2(1) within different error bars, corresponding to 
different observables measured during simulations. The discrepancy between 
our estimates an those of Refs. [34,26] may be caused by the fact, that the 
latter have defined scaling exponents for the out-of-equilibrium short-time dy- 
namics. 



The Swendsen-Wang and Wolff algorithms give rise to the dynamics of spin 
clusters, which differs from the local one. Even for the pure 3d Ising model 
there is no field theory describing such dynamics. However, there exist esti- 
mates, relating dynamical critical exponent of the cluster algorithms to the 
static exponents. Besides Eqs. (20), (21), the following inequality has been 
proven for the energy-like integrated and exponential autocorrelation time 
critical exponents of Swendsen-Wang algorithm [22]: 



^SW ^SW ^ ^, /,, fr,o\ 



Eqs. (20), (21), (23) hold for the pure Ising model. In particular, Eq. (23) leads 
to the conclusion, that systems with a positive specific heat exponent a must 
display critical slowing down. In absence of such inequality for the diluted 
system, our results for the critical dynamics of RIM for cluster algorithms 
prove that the critical slowing down is present in diluted systems as well. 
However, a striking feature of our estimates (second and third lines of table 6) 
is that they suggest that dilution leads to decrease of the dynamical critical 
exponent for the cluster algorithms. This phenomena is quite opposite to the 
local dynamics, where dilution enhances critical slowing down. The values of 
the exponents describing relaxation of different observables differ numerically, 
being however close to each other. Nevertheless, on this stage it is impossible 
to exclude that the difference is not only due to the crossover phenomena [49] . 
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